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Richardson approach provides an exact solution of the pairing Hamiltonian. This Hamiltonian 
is characterized by the electron-hole pairing symmetry, which is however hidden in Richardson 
equations. By analyzing this symmetry and using an additional conjecture, fulfilled in solvable 
limits, we suggest a simple expression of the ground state energy for an equally-spaced energy-level 



model, which is applicable along the whole crossover from the superconducting state to the pairing 
fluctuation regime. Solving Richardson equations numerically, we demonstrate a good accuracy of 
£f} • our expression. 
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vL| ! I. INTRODUCTION 

Or 

CZ3 ' As shown by Richardson, the "reduced" Hamiltonian of the Bardeen-Cooper-Schrieffer (BCS) theory of super- 

conductivity is exactly solvable*. The Richardson approach is based on the canonical ensemble, i.e., the number of 
particles is fixed, while the BCS theory corresponds to the grand canonical description, which becomes accurate in the 

i-^ ' large-sample limit. In addition, the BCS theory provides only a mean-field approximation, which nevertheless turns 

out to be exacl^ in the large-sample limit due to peculiarities of the "reduced" BCS potential. However, the grand 

^j • canonical BCS theory is definitely not applicable for small-sized systems accommodating few pairs only. Significant 

improvements, nevertheless, are possible if one incorporates canonical approach into the BCS model (particle-number 

projected BCS), but still a rather heavy numerics is needed to proceed with computations^. For the applicability of 

^^ . the canonical ensemble to the theory of superconductivity see Ref. 4 . 

In contrast to the BCS theory, the Richardson method yields an exact solution to the many-body problem involving 
BCS pairing Hamiltonian. Within this approach, the energy of the system of N correlated pairs is expressed through 
the sum of N energy-like quantities, which satisfy the system of N coupled nonlinear equations, called Richardson 
equations. It is remarkable that they are also derivable through the algebraic Bethe-ansatz approach^, so that the 

t-H | Richardson equations can be considered as one of the examples of Bethe-ansatz equations. However, it turns out 
^ \ that solving the Richardson equations is a formidable task. Up to now, they have been evaluated explicitly only in 

S^ few special cases. In particular, no analytical solution exists for the crossover between the superconducting state and 

the pairing fluctuation regime, which is relevant for small-sized systems even for zero temperature. In this situation, 
no small parameter seems to exist, which could be used to construct some expansion. For this case, the Richardson 
equations are tackled numerically, when studying small systems described by pairing Hamiltonian, among which are 
nanosized superconducting grains®, nuclei 7 , ultracold atoms^, bubbles in liquid helium 9 (for studies of correlation 
functions, see Refj^). 

The aim of the present paper is to apply symmetry arguments in order to provide analytical results for the crossover 
region. It is actually well known that symmetry considerations can be very helpful in situations, when brute-force 
methods are not so efficient. We show that the BCS pairing Hamiltonian has an electron-hole (e-h) symmetry and 
we then try to reveal the impact of this symmetry for the ground state energy for the arbitrary number of pairs and 
interaction constant. The case of the equally-spaced model is considered, when the energy levels of noninteracting 



particles are distributed equidistantly. The impact of e-h symmetry can be most readily revealed in this particular 
case. 

Although the e-h symmetry is encoded in the pairing Hamiltonian, it does not show up explicitly in Richardson 
equations. From the symmetry arguments, we suggest a simple explicit formula for the ground state energy, which 
constitutes a main result of this paper. Within our approach, we use a conjecture on the TV-dependence of the 
dominant contribution to the ground state energy (when all other parameters are fixed). This hypothesis is justified 
by the fact that a guessed simple dependence on N emerges in very different analytically-solvable limits. Under this 
assumption, we actually reduce the problem to the resolution of a single Richardson equation, which is far simpler 
than the solution of the full set of equations. The resulting expression for the ground state energy is also in a perfect 
agreement with known results in these limits, although we do not solve the whole system of Richardson equations. 
In order to address the accuracy of our formula in the crossover regime, we solve numerically the full system of these 
equations for N pairs changing from N — 1 to 50 and compare the results. We found a very good accuracy for systems 
with small number of pairs, which are of particular interest due to limitations of BCS treatment in this case. For 
systems with larger number of pairs, the maximum error, which corresponds to the regime of intermediate strength 
of coupling, grows. Nevertheless, our expression remains applicable in this case too. Our approach can be useful for 
other types of Bethe-ansatz equations, especially for those, which correspond to Gaudin-like models. 

Though BCS Hamiltonian is rather simple, its applicability to small diffusive superconducting grains has been 
proved, see, e.g., Refs.— ~— . The crucial condition is that the Thouless energy, which gives the inverse time to diffuse 
across the grain 14 , must be much larger than the average energy level spacingii. Also important is to have the 
Thouless energy much larger than the superconducting gap (see Appendix B of Ref. 12 ). It is less obvious what should 
be a relation between the energy level spacing and the superconducting gap. In particular, it was argued in Refi 1 ^ 
(see also Ref.—) that the BCS Hamiltonian probably can be considered as a toy model only, when the gap becomes 
much smaller than the level spacing. 

II. MODEL 

We consider a system of fermions of two sorts, for instance, with spins up and down. Particles attract each other 
through the BCS "reduced" potential, coupling only fermions of different sorts and with zero total momenta as 

V = -V^ a k't -k'4 a -HOkf. (1) 

k,k' 

The total Hamiltonian is H = Hq + V, where 

H = ^ £k f OkfOkt + a k|. a kU • ( 2 ) 

k 

The summation in the right-hand side (RHS) of Eq. (fTJ runs only over the states with kinetic energies £k and £k' 
located in the energy band between sp and £f + ^ (Debye window). These energies are distributed equidistantly 
(equally-spaced model), so that the difference between two nearest values of e^ is 1/p. The density of states p increases 
with the system volume, while the interaction constant V decreases, so that the dimensionless interaction constant 
v = pV in superconductors is finite and, in the large-sample limit, it can be treated as a material characteristics. In 
the BCS theory, the energy interval between £f and £p + Q is assumed to be always half-filled, while ft/2 is the 
Debye frequency. Thus, the total number of available states with up or down spins in the Debye window is Nq = pQ, 
while N = Aq/2; eu runs from £p a to Ef + (Nq — l)/p taking A$i values in total. 

While it is supposed that the half-filling configuration only is physically meaningful, one can consider other fillings 
and at least treat the problem from the purely mathematical perspective. Note however that it was argued in Refj 1 ^ 



that such a model might be relevant to some semiconductors (see also RefJ£). Introduction of the extra degree 
of freedom, which is a filling, is an important ingredient of our analysis. By considering the energy of the system 
formally as a function of N (with all other input parameters fixed) , we are going to obtain a valuable information on 
the half-filling situation. Numerical results will be presented for N = Nq/2 only. For the sake of simplicity, we will 
focus on even values of Nq. 

In the present paper, we concentrate on the equally-spaced model, which provides a simplest but physically mean- 
ingful distribution of energy levels. Therefore, this model is the most attractive starting point to study the impact 
of e-h symmetry on the solutions of Richardson model. Superconducting correlations were also studied in finite-size 
system with random spacings of levels, distributed in accordance with the gaussian orthogonal ensemble, using mean- 
field 1 ' and exact Richardson approaches 18 . Such a statistics is typical for small metallic grains 19 . It was shown 18 
that the crossover between the superconducting state and the fluctuation-dominated regime is smooth, similarly to 
the case of the equally-spaced model. 

The Hamiltonian, defined in Eqs. (JT|) and ([2]), is exactly solvable^. The energy of N pairs is given by the sum of N 
rapidities Rj (j = 1,..., N) as En = Ylj Rj- The Richardson equation for each rapidity Rj reads 

where the summation in the first term is performed over £k located in the Debye window. Note that the dependence 
of En on N thus enters through the number of equations. 

III. THE ELECTRON-HOLE SYMMETRY 

Now we discuss the internal electron-hole symmetry contained in the Hamiltonian, from which we are going to deduce 
an information on the solutions of Richardson equations. Let us introduce creation operators for holes as b^ = aicf- 
and frj. = akj.. By using commutation relations for fermionic operators, it is easy to rewrite the Hamiltonian in terms 
of holes as 



H = -VNn+ 2 Y, £ *-Y,( £ *- V l( b U bk t + b tM) 

k k 

~ V E^-k^-H&kf (4) 



k,k' 



The first two terms of the RHS of Eq. (j4]) are numbers. They give the potential energy and the kinetic energy 
of the Debye window completely filled by electron pairs, respectively. The fourth term coincides precisely with the 
interaction potential in terms of electrons, given by Eq. (fTJ). To analyze the third term, we introduce £ k , defined as 
£ k = Ef„ + (Nq — l)/p — £k, which takes values 0, 1/p, 2/p, ..., (Nq — l)/p, so that £ k runs over all the states starting 
from the top of the Debye window towards its bottom, i.e., in the inverse order. Then, — (ek — V) can be represented 
as £ k + (V — sp — (Nq — l)//o). A similar term in the Hamiltonian for electrons, given by Eq. ([2]), contains a factor 
£k + £ f , where £k takes values 0, 1/p, 2/p, ..., (Nq — l)//», so that it runs over all states starting from the bottom 
to the top of the same energy band. 

Thus, due to peculiarities of the interaction potential, there exists a symmetry between electron and hole pairs in 
the Hamiltonian. Moreover, due to the equally-spaced distribution of energy levels in the Debye window, the ground 
state energy of N electron pairs can be explicitly expressed through the energy of Nq — N electron pairs, with ef 
changed into (V — £f — (Nq — l)/p)- We therefore treat this energy En(ef ) as a function of Ef , which is an arbitrary 
nonzero number, and a discrete variable N, which runs over the set 1,2, ...,Nq. Since the bare kinetic energy of a 



filled Debye window is given by the sum of terms of the arithmetic progression, 2Nqe Fo + Nq(Nq — l)/p, we arrive 
at the identity 

Nn-V 



En(£Fo) — ^Nn-N [V - e Fo 



P 



-VN n + 2N n s Fa + Nn{Nn - 1) . (5) 

P 

Note that in the case of more complex types of energy-level distributions in the conduction band, one has also to 

change the distribution, when switching from electron pairs to hole pairs, which corresponds to counting levels from 

the top instead of counting them from the bottom of the Debye window. This makes the situation more subtle 

compared to the equally-spaced model. Nevertheless, the duality between the electrons and holes still exists, since 

this feature is a direct consequence of BCS interaction potential, so that Eq. (jj) stays valid. 

Next, we split En{£ Fo ) into the additive contribution 2Ne Fo , which simply corresponds to the shift of all rapidities, 

and E N , the latter being independent of e Fo 

E N (e Fo )=2Ne Fo +E' N . (6) 

By substituting Eq. ([5]) to ([5]), we arrive at the functional equation 

iVn-1 



E N = Eno-n + (Nn - 2N) \V - -^— J , (7) 

which is still exact. This remarkable condition will enable us to relate the condensation energy of N pairs to the 
condensation energy of Nq — N pairs. Note that it is automatically fulfilled for the half-filling, N = Nq/2. 

We would like to stress that, although the e-h symmetry can be rather easily extracted from the Hamiltonian, 
it is not obvious that the solution of N equations is related to that of Nq — N equations through such a simple 
and universal relation as Eq. ([7|). For the moment, we do not know how this relation can be derived directly from 
Richardson equations. 

Since E N is a function of the discrete variable N , which runs over Nq values, it can always be represented as 
a polynomial of N of power Nq. Alternatively, instead of expanding in elementary monomials N n , one may use 
Pochhammer symbols defined as 

(N) n = N(N-l)...(N-n + l), (8) 

while (N)o = 1; so that (N) n may be treated as a polynomial of N of power n. Then, 

E N = J2 R n(N) n , (9) 

71=1 

where e„ is a set of unknown numbers. 

Actually, E N can be also split into the condensation energy E^ on and the contribution coming from the bare 
kinetic energy. The latter is given universally by N(N — l)/p, which can be obviously described by the second term 
in the RHS of Eq. ©. 

Up to now, all the results were exact. At this step, we make a conjecture that a dominant contribution to E N is 
due to the first two terms in the sum of the RHS of Eq. (j9|) that is 

E N ~e x N + eiN{N-l). (10) 

This assumption is fully reasonable, since such a form of E N does emerge in three important limits, which are solvable 
analytically. Let's discuss the condensation energy in these three limits, since the contribution from the kinetic energy, 



N(N — 1")/ p, is always in agreement with Eq. flTOl) . as discussed above. The first limit is a regime of the very weak 
coupling realized for finite Nq (v <C l/ln(JVn)), for which all the rapidities are located in real axis and approach the 
energy levels of noninteracting electrons. In this case, E^ on — —VN, so that Eq. (jlOl) is satisfied. Another limit is 
the strong-coupling regime (v 3> 1), when all the rapidities are located far away from the line of one-electron levels 
in the complex plane. In this case, Ej^ on contains 1 terms proportional to N and N(N — 1). At last, there is a limit 
of infinite N at finite nonzero u; here also E^° n has a similar form 2 ^. These three limits are quite different from 
each othe r 18 ! 21 ' 22 . It is actually the reason why we conjecture that the structure of the solution, given by Eq. (|10j) . 
must remain robust in the intermediate region, which is characterized by finite Nq and arbitrary v. 

Now we consider e± and e^ as unknown numbers and substitute Eq. (fT0|) into Eq. ([7|). We then equate coefficients 
of (iV)o, (N)i and (N)2 in both sides of this equation and obtain a system of three linear equations for e\ and e-i- 
These three equations turn out to be dependent, so they yield only a single condition as 

gi , f Nn-l \ 1 f . 

e2 = -^-l + {— V )^-l- (U) 

At this stage, we are left with only one unknown number, e±. It can be easily determined by considering a one-pair 
problem (as a "boundary condition" in the space of discrete N). In this case, E 1 = ei, while e\ is given by the 
solution of the single Richardson equation 

Nn-l . 

under the condition e\ < 0, which ensures that we select the lowest-energy solution; then, —e\ is a binding energy of 
a single pair. In the general case, e\ has to be determined numerically, while exact analytical results are available in 
certain limits, as discussed below. Note that Eq. (fT2|) can be rewritten in terms of r-functions. 
The expression of E N is obtained by substituting Eq. (fTTj) to Eq. (fTOj) as 

N(N-l) ( N-l\ N-l 

where the first term comes from the bare kinetic energy, while two others give the condensation energy. It is remarkable 
that our method reproduces the first contribution automatically in the exact form. The condensation energy per pair 
then reads 

Eq. (UJ| supplemented by Eq. (fTi2)) is the main result of this paper. 



IV. DISCUSSION 



Let us now consider the limits when Eq. (|12|) can be solved analytically, in order to see whether Eq. (jl4|) gives 
reasonable results. 

We start with the limit of the very weak coupling, in which the solution of Eq. (|T2j) approaches the lowest level so 
closely that the mutual separation becomes much smaller than 1/ ' p and, moreover, contributions from other levels can 
be neglected. Then, we obtain a simple solution e\ = —V. By estimating dropped contributions due to other levels, 
we obtain a criterion of applicability of this result as v <C l/ln(7Vn). In this case, E^ on ' /N reduces to —V, as it 
must be, due to the cancellation in the RHS of Eq. (fl4|) . Note that the condensation energy per pair in this limit is 
independent of filling. 



Next, we consider an opposite limit, when the separation between the single rapidity and the lowest level is much 
larger than 1/p and the number of levels Nq is also large. This enables us to replace the sum in the RHS of Eq. (TT2")) 
by the integral, which gives 

ei c-2n -^tlM (15) 

l-exp(-2/u) v ; 

The condition of applicability of this result is thus twofold: Nq ;§> 1 and Nq exp(— 2/v)/{\ — exp(— 2/v)) ^> 1. The 
infinite-sample limit, when v is fixed and finite, thus always satisfies these criteria. It is also easy to see that | e\ \ 
given by Eq. (|T2|) is much larger than V in the large-sample limit, Nq — > oo, so that V ~| e\ \ /Nq. This means 
that in this limit V can be neglected in Eq. (fT4|) . Then, 



Acond) , M ^_ fr>_ N\ exp(-2/u) 

p ) l-exp(-2/«)' 



E^/N ~ -2 ( n - - ) i ^Zi.r.L : (16) 



which, for the half-filling, coincides with the BCS expression for the condensation energy per N (and with Richardson 
large- N result for the same filling 2 ). For arbitrary filling, it also coincides with the results of both the mean- field 
treatment and of the Richardson approach 2 ^. 

We would like to stress that the obtained results are highly nontrivial. Indeed, the condensation energy given by 
Eq. (|14p consists of two terms. In the limit of a very weak coupling, both of them are of the same order, while their 
combination gives the exact result within all numerical prefactors. In contrast, in the large-sample limit Nq — > oo, 
when interaction constant v is finite and nonzero, one of the terms becomes of the order of 1/Nq compared to 
another one, so that it can be dropped as an underextensive contribution, while the remaining term again gives a 
correct result within all numerical prefactors. Such a very subtle interplay between the two terms indicates that the 
suggested formula for the condensation energy must remain accurate not only in the considered limits. 

Actually, by manipulating the single-pair binding energy, which must be determined numerically, we circumvent 
the problem of inapplicability of large- N approaches to the normal state, pointed out by Richardson in Ref3. This 
difficulty can be traced back to the non-analytic dependence of condensation energy on V in this limit, while in the 
limit of a very weak interaction it is simply proportional to V. Moreover, in the first case, this energy is an extensive 
quantity, while in the second one, it becomes intensive. It is therefore challenging to unify both regimes within a 
single self-consistent formalism. 

In order to explore the applicability range of the obtained expression of the condensation energy, we perform a 
systematic numerical solution of the full set of Richardson equations, as given by Eq. ((3|), for various values of N from 
1 to 50 and v at the half-filling. Then, we compare the obtained numerical results with the prediction of Eq. (Til)) 
(where e\ is obtained from the numerical solution of Eq. (|12j) ). We also calculate the condensation energy by using 
the standard grand-canonical BCS theory. The only difference with the common version of this theory is the fact 
that we do not replace sums by integrals when solving the gap equation and also when calculating the condensation 
energy itself, which is of importance for systems with relatively small number of pairs, since these replacements are 
responsible for additional inaccuracies. 

In the general case, the Richardson equations can be numerically solved by Newton-Raphson method with a good 
initial guess. An exact solution Rj = 2e^ Q + 2j/p (n — 0, 1, .., N — 1) is a solution for v — 0. Therefore, we start with 
such initial values and then find solutions with increasing v. In order to avoid the singularity, new variables A+,A_ 
are introduced^. When v is close to the critical v c , the Newton-Raphson result does not converge to the solution if 
it starts from the other side of the singularity 24 . Therefore, an extrapolation step is taken for the new v (close to v c ), 
as proposed in Rcf. 24 . 

The results are presented in Fig. 1 for three particular values of N. The condensation energies are measured in 
terms of 2/ p. Fig. l(a - c) give the dependence of the condensation energy per pair, as a function of v, for N = 5, 25, 
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Figure 1: (Color online) The condensation energy per pair as a function of the interaction constant, for TV pairs: TV = 5 (a), 
TV = 25 (b), and TV = 50 (c). Red (gray) solid lines show the prediction according to our analytical formula (|14| . blue (light 
gray) dashed lines represent the results of the numerical solution of the full system of Richardson equations, black dotted lines 
correspond to the grand-canonical BCS result. 



and 50, respectively. Solid curves yield our prediction, dashed curves represent results of the numerical solution of the 
Richardson equations, and dotted lines are the grand-canonical BCS results. We see that there is generally a good 
agreement between the numerical results and those obtained from our formula. The similar agreement has been found 
for other values of TV. Therefore, our conjecture is justified. In contrast, BCS results become accurate in the large- TV 
limit only: as it is known, there is a range of small v for any TV, when the grand-canonical BCS theory is qualitatively 
incorrect, since it predicts a disappearance of superconducting correlations (see Fig. 1). In this fluctuation-dominated 
regime, our approach works well. We found that the largest relative errors for the three cases illustrated in Fig. 1 are 
3, 15, and 23 percent, respectively. Thus, our approach is more efficient for small TV, this case being more interesting 
due to limitations of BCS theory for systems with small number of pairs. We also would like to note that the accuracy 
can, in principle, be improved by considering more terms in Eq. ^ and using two-pairs "boundary condition" (and 
possibly other configurations with even more pairs), although such a procedure would lack simple analytical results, 
in contrast to the present approach. 
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Within our approach, a single-pair binding energy —e\ plays a very important role, although we deal with the 
many-pair system. It provides an energy scale, which is alternative to the superconducting gap A. Interestingly, the 
existence of an additional scale for finite systems was revealed some time ago in Refill, where it was shown that 
BCS results for the ground state energy become inadequate already for level spacings 1/p « 2A 2 /$7, which are much 
smaller than A at « < 1. Let us point out that, as easily seen by a direct comparison, in the large-sample limit, 
2A 2 /f2 is nothing but — e\ when v is small. It is very unlikely that such a coincidence is accidental. Therefore, we 
believe that the additional energy scale found in Rcf. 21 is connected to the single-pair binding energy. Note that it 
was recently argued 2 ^ that standard BCS results for the condensation energy in the thermodynamical limit can be 
also interpreted in a simple way through the single-pair binding energy and not through A, as usual. In particular, 
in this limit, — e\ and A have similar, but different dependencies on v, since — e\ ~ exp(— 2/v), while A ~ exp(— 1/v) 
at v <C 1, as discussed in detail in Ref^. 

It is perspective to extend our analysis to excited states. We also think that, probably, similar ideas can be applied 
to other types of energy-level distributions, in particular, to those, which are characterized by the inversion symmetry 
around the middle point of the band. Furthermore, the analysis of the electron-hole symmetries can be helpful in 
treatments of other Bethe-ansatz equations, among which Richardson equations are just one of the examples. 

Notice that the derived exact relation for the condensation energy, given by Eq. (0 , can be used as a test to check 
the accuracy of the numerically-found solutions of Richardson equations. 

V. CONCLUSIONS 

We suggested a simple expression for the ground state energy of the pairing Hamiltonian for the case of the equally- 
spaced model along the crossover from the superconducting regime to the pairing fluctuation regime. This expression 
is derived from the peculiar electron-hole symmetry of the pairing Hamiltonian and relies on the conjecture, which 
enables us to reduce the task to the one-pair problem. The electron-hole symmetry is encoded in the Hamiltonian, 
but it is hidden in Richardson equations, which provide an exact many-body solution of the problem. The obtained 
expression of the ground state energy depends on the binding energy of a single pair, which, in the general case, must 
be determined numerically by solving a single Richardson equation. The latter problem is much simpler than the 
solution of the full set of Richardson equations. This quantity also provides an additional energy scale associated with 
superconducting correlations, which seems to be connected with the energy scale revealed in Refill. 

The comparison with the results of the full numerical resolution of Richardson equations demonstrated a generally 
good accuracy of the suggested formula, while the usual grand-canonical BCS approach fails even qualitatively in the 
fluctuation-dominated regime. The accuracy is better in the case of the system with small number of pairs, which is 
of particular interest due to limitations of BCS method in this situation. 
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